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ABSTRACT 



> 

' Molecular-line observations of star- forming cloud cores indicate that they are not the flattened 

. structures traditionally considered by theory. Rather, they are elongated, perhaps in the direction 

of their internal magnetic field. We are thus motivated to consider the structure and evolution 
of axisymmetric, magnetized clouds that start from a variety of initial states, both flattened 
' (oblate) and elongated (prolate). In this first contribution, the clouds are of fixed mass, and are 

surrounded by a fictitious medium of zero density and finite pressure. We devise a new technique, 
O \ dubbed the g-method, that allows us to construct magnetostatic equilibria of any specified shape. 

The mass loading of the field lines then follows from the self-consistent model solution, just the 
^ ' reverse of the standard procedure. We find, in agreement with previous authors, that the field 

lines in oblate clouds bend inward. However, those in prolate clouds bow outward, confining the 
structures through magnetic tension. 
, We next follow the quasi-static evolution of these clouds via ambipolar diffusion. An oblate 

' cloud either relaxes to a magnetically force-free sphere or, if sufficiently massive, flattens along 

its polar axis as its central density runs away. A prolate cloud always relaxes to a sphere of 
modest central density. We finally consider the evolution of an initially spherical cloud subject 
to the tidal gravity of neighboring bodies. Although the structure constricts equatorially, it 
also shortens along the pole, so that it ultimately flattens on the way to collapse. In summary, 
none of our initial states can evolve to the point of collapse while maintaining an elongated 
shape. We speculate that this situation will change once we allow the cloud to gain mass from 
its environment. 

Subject headings: ISM: clouds — ISM: structure — ISM: magnetic fields — MHD — stars: formation 

1. INTRODUCTION ture of the dense, molecular cloud cores that form 

low-mass stars (see, e.g., Myers 1999). A typi- 
After two decades of intense observational ef- ^^1 core, with a mass of several solar masses and 

fort, astronomers have learned much about the na- ^ diameter of order 0.1 pc, has a central density 

; somewhat in excess of 10^ cni^^ and a kinetic tem- 
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still only margmally feasible at such densities, sug- 
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gest magnetic field strengths no larger than 30 /iG 
(Crutcher 1999). These figures refer specifically 
to dense cores in relatively sparse environments 
like Taurus- Auriga. Those within regions such as 
Orion, that form high-mass stars, appear to be 
more massive and warmer, although the measure- 
ments here are generally less precise because of the 
greater distances involved (e.g., Jijina, Myers, & 
Adams 1999). 

To see how a dense core evolves to form a star, 
one must first have a physical understanding of its 
structure. A core's internal velocity dispersion, 
as gauged from the linewidths of tracer molecules 
like NH3 and CS, depends on its mass. In low- 
mass cores, the dispersion is typically 0.2-0.4 km 
s~^, significantly lower than the background gas. 
This range in dispersion matches that expected 
from the virial theorem for objects of the appro- 
priate mass and size. The line profiles, moreover, 
show only slight nonthermal broadening (Myers 
& Benson 1983). Finally, the internal magnetic 
field strength is also consistent with the virial ex- 
pectation, and with approximate equipartition of 
thermal and magnetic energies (Myers & Good- 
man 1988). These facts together imply that the 
typical low-mass core is in dynamical equilibrium, 
supported against self-gravity by a combination of 
ordinary gas pressure and the Lorentz force asso- 
ciated with a largely static magnetic field (McKee 
et al. 1993). 

Theorists have long recognized that such a 
magnetostatic structure evolves quasi-statically 
through ambipolar diffusion, i.e., the relative drift 
of neutrals and ions (Mestel & Spitzer 1956). As 
the configuration slips through the ambient mag- 
netic field, its central density gradually increases 
to the point where rapid, protostellar collapse 
begins (Nakano 1979). While this basic picture 
continues to provide a framework for our con- 
ception of dense core evolution, its quantitative 
implementation is not without problems. One im- 
portant issue, and the focus of this paper, concerns 
the three-dimensional shape of the objects. 

Radio maps in a number of tracer lines show 
dense cores to be distinctly nonspherical, whether 
they are in loose associations (Benson & My- 
ers 1989) or massive, turbulent complexes (Harju, 
Walmsley, & Wouterloot 1993), and whether they 
contain stars or not (Jijina et al. 1999). The ob- 
served distribution of projected axial ratios is well 



matched if the cores arc intrinsically prolate struc- 
tures, with a random orientation of their long axes 
in the plane of the sky (Myers et al. 1991; Ry- 
den 1996). Theoretical models of magnetostatic 
clouds, on the other hand, generally predict that 
they are oblate, flattened along the direction of the 
background magnetic field (Mouschovias 1976b; 
Tomisaka, Ikeuchi, & Nakamura 1988a). Under 
ambipolar diffusion, the central region flattens 
even more as the central density climbs (Tomisaka, 
Ikeuchi, & Nakamura 1990). 

A useful step toward reconciling theory and ob- 
servation is to broaden the range of models and 
consider cores with a variety of intrinsic shapes, 
including prolate configurations. This is the first 
task of the present study. The essential reason 
why previous authors obtained flattened struc- 
tures is not difficult to see. In solving the equa- 
tions for force balance, one customarily specifies 
the amount of mass loading each magnetic fitix 
tube (Mouschovias 1976a). Lacking any obser- 
vational indication of this mass-to-flux distribu- 
tion, theorists have taken it from geometrically 
simple reference states, such as uniform spheres 
or cylinders threaded by a spatially constant field. 
Such configurations, however, are clearly not in 
force balance. If allowed to evolve, they imme- 
diately collapse along the field until dynamical 
equilibrium is restored (see, e.g. Fig. 2 of Fiedler 
& Mouschovias 1993). The magnetostatic struc- 
tures resulting from such collapse are inevitably 
flattened, and become more so during the subse- 
quent phase of quasi-static settling via ambipolar 
diffusion. This qualitative behavior is unaltered 
by the presence of isotropic, subsonic turbulence 
(sec, e.g., the schematic treatment of Lizano & Shu 
1989). 

To remedy this sitiiation, we break from tradi- 
tion and use the core shape itself as input, rather 
than the unknown mass-to-flux distribution. That 
is, we solve the magnetostatic equations subject to 
the constraint that the gas pressure fall to some 
uniform, ambient value along a specified bound- 
ary. The mass-to-flux distribution then results 
from the self-consistent numerical solution. We 
use this technique, dubbed the q-method, to con- 
struct a variety of equilibria, with shapes rang- 
ing from oblate to prolate. We further distin- 
guish two classes of states, those which are spa- 
tially isolated, and those subject to the tidal grav- 
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itational field of nearby cores. Study of the first 
class allows ready comparison with previous au- 
thors, while the tidal boundary condition approx- 
imates the observed situation within filamentary 
clouds, where neighboring cores are typically sep- 
arated by only 0.5-1.0 pc (e.g., L1495 in Taurus; 
Onishi et al. 1996). 

Our next task is to follow the quasi-static evo- 
lution of selected initial states through ambipolar 
diffusion. In this paper, we assume each core to 
be of fixed total mass. That is, we take the am- 
bient pressure to arise from a fictitious medium of 
negligible mass density and correspondingly high 
temperature. Our numerical results first confirm 
that isolated, oblate structures undergo significant 
central flattening as their density rises steeply. Iso- 
lated, prolate structures experience no runaway 
increase in density, but instead relax to spheres. 
Finally, an initially spherical, but tidally stressed 
configuration does undergo central condensation, 
but concurrently shrinks along the polar axis. In 
summary, none of our states both evolves to a high 
central density and maintains an elongated shape, 
as the observations seem to demand. 

Section 2 below details our numerical method 
for constructing magnetostatic equilibria with 
specified shapes. In Section 3, we describe the 
physical properties of these equilibria. Section 4 
then presents numerical results for the quasi-static 
evolution of a few representative models. Finally, 
Section 5 assesses our results in light of previous 
observational and theoretical work. We conclude 
by emphasizing the potential value in relaxing the 
assumption of fixed core mass. In a second paper, 
we shall explore cloud evolution under these more 
general circumstances. 

2. FORMULATION OF THE PROBLEM 
2.1. Governing Equations of Equilibrium 

The theory of self-gravitating, magnetized 
cloud equilibria, first discussed by Mestel (1965) 
and Strittmatter (1966), has been extended 
and reformulated by many authors notably 
Mouschovias (1976a), Nakano (1979, 1984), and 
Tomisaka et al. (1988a). Here, we briefly review 
the basic equations, before passing to a discussion 
of the features unique to our method. 

The governing equations are the condition of 
force balance, Poisson's equation, and Ampere's 



law: 

-VP-/)V^ + -jxB = 0, (1) 

c 

VV = 47rGp, (2) 

VxB = — j, (3) 

c 

where the chosen symbols are standard. Our 
adopted equation of state is that of an isother- 
mal gas, P — a^p, where a is the constant sound 
speed. The ancillary relation 

B = V X A, (4) 

where A is the magnetic vector potential, ensures 
that V-B = is satisfled identically. 

We adopt a cylindrical coordinate system 
(r, 4>, z) whose origin lies at the cloud center; see 
Figure 1. We assume axial symmetry about the 
2;-axis, and reflection symmetry about the plane 
z = 0. Hence, only the quadrant r > 0, z > 
needs to be considered explicitly. We assume that 
B is poloidal, and hence A toroidal. We thus take 
A = A{r, z) 60. Alternatively, B = — r^^e^ x V$, 
where $ is the scalar function defined by 

^{r,z) = rA{r, z). 

The quantity $ is invariant along each magnetic 
surface: B-V<I> = 0. Since it is proportional to the 
usual magnetic fiux i.e. $ = ^b/Stt, we refer 
to loosely as the "magnetic flux," and use it to 
label field lines. 

It is convenient to resolve the forces in equa- 
tion (1) both along and across field lines. Follow- 
ing Dungey (1953) and Mouschovias (1976a), we 
introduce the scalar fmiction 

g ee; P exp(V'/a^) = o^p exp(i/;/a^). (5) 

Equation (1) may be rewritten in terms of $ and 
q as 

j V$ = cr e^i-ip/a^) Vq, (6) 

where we have noted that j = j e^. While q de- 
pends on both r and z, its real utility comes from 
the fact that it is a function of $ alone. That is, 
it follows from equation (6) and B-V<I> = that 
B-Vg = 0. This property represents force balance 
in the direction parallel to B. Since q = q{^), 
we may use equation (6) to express force balance 
perpendicular to the field as 

i:exp(^/a^) = (7) 
cr el's 
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Using equations (4) and (7), equation (3) becomes 



d_ 
dr 



ld_ 
r dr 



{tA) 



—A-KT exp{—'tp/a'^ 







(interior) (8) 
(exterior), 



where the terms "interior" and "exterior" refer to 
the regions inside and outside the cloud, respec- 
tively. As in previous studies, the cloud exterior 
is assumed to be composed of a hot and tenuous 
gas that has no dynamical effect except to exert 
a finite pressure, Pq, on the cloud boundary. The 
exterior region is itself bounded by a cylinder of 
radius R and half- height Z (sec Figure 1). It is on 
this cylindrical surface that we apply our bound- 
ary conditions (see §2.2 below). 

Equation (8) is to be solved throughout the en- 
tire volume, simultaneously with Poisson's equa- 
tion (2), 



ld_ 
r dr 

AttG 



dr J dz"^ 
qi^) exp(-V'/a2) 




(interior) (9) 
(exterior). 



where we have made use of equation (5) on the 
right-hand side of the interior equation. Solving 
the equations in this form guarantees that the rel- 
evant continuity conditions for the gravitational 
and magnetic fields across the cloud surface are 
satisfied (Mouschovias 1976a). 

2.2. Boundary Conditions 

The boundary conditions are applied along the 
cylindrical surface bounding the exterior zero- 
density region. The boundary condition on A is 
that the magnetic field approach a uniform back- 
ground value, B = Boo, at the boundary. That 
is, 



A = 



2 

BooR 



at z = Z, 
ad, r = R. 



(10) 



We consider two types of boundary condition on 
the potential tjj, reflecting different assumptions 
regarding the external medium. 



(i) Point-Mass, or "Isolated" Boundary Condi- 
tion. This boundary condition has been adopted 
in most previous studies. It asserts that tp in the 
exterior region far from the cloud is indistinguish- 
able from that of a point particle of the same mass: 



V'(r, z) = 



-GM 



(^2 _^ ^2)1/2 



at r = R, z = Z, (11) 



where M is the mass of the cloud. In practice, 
we have used R = 2Rq and Z = 2Zq when ap- 
plying this boundary condition. Here Rq and Z^ 
are, respectively, the radial and vertical extent of 
the cloud. Calculations employing larger cylinders 
give essentially identical results. 

(ii) Periodic, or "Tidal", Boundary Condition. 
This boundary condition dictates that the gravi- 
tational force vanish along some surface at a finite 
distance from the cloud, as is appropriate when 
surrounding matter is present. We take these sur- 
faces to be the top and bottom of our cylindrical 
outer boundary (e.g. Fiedler & Mouschovias 1992). 
The conditions on the potential then become 



di}_ 
dz 



at 2; = Z, 
constant at r = R. 



(12) 



Note that the first condition holds if the cloud of 

interest is one of an infinite chain, with a spacing 
of 2Z. The constancy of the potential at r = i? 
simulates the presence of a larger background fil- 
ament.^ 

Finally, the assumed refiection symmetry about 
= implies that 



dz dz 



at z = Q, 



while the vanishing of the radial gravitational force 
and the lack of magnetic sources along the 2;-axis 
imply 

^ = 0, A(0, z) = Q at r = 0. 
dr 



^Our periodic boundary condition differs from that used by 
Lizano & Shu (1989), in that the latter authors took the 
gravitational potential to be constant on a specified locus 
corresponding to a chain of point masses. Moreover, they 
forced B to equal a uniform background value on that sur- 
iaice. 
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2.3. Construction of Initial States: 

g Method 



The 



Consider a single magnetic flux tube penetrat- 
ing the cloud, containing mass 6m and flux 6^. 
The former is given by (Mouschovias 1976a) 



(5m($) 



Jo Jr(<S> 



r(3>+5$,2) 



2Trrp{r, z)dr, 



where Zci(r) describes the boundary shape (Figure 
1). Changing variables from r to $, we have dr = 
d^{dr/d^), so that a trivial integration over $ 
yields the cloud's mass-to-flux distribution 



dm 



X exp 



/ r($,2;) 
Jo 



dz, 



(13) 



where we have substituted for p from equation (5). 
Here m($) is the total mass contained within a 
given (axisymmetric) surface of constant magnetic 
flux 

dm($) 



;■»,($) 



d^ 



d^. 



(14) 



Denoting the total fltix enclosed in the cloud by 
$0 = *(-Ro), the total mass is then M = m($o)- 

The customary procedure for obtaining magne- 
tostatic equilibria is to solve equation (13) for q, 
obtaining 



9($) = 



X exp 




(15) 



This expression is then used in the right-hand sides 

of equations (8) and (9). Note that all quantities 
on the right-hand side of equation (15) are deter- 
mined iteratively in a numerical scheme, except 
dm/d^, which is specifled a priori (§1). 

As noted in §1, we invert the customary pro- 
cedure when constructing our initial states. We 
do not specify a priori dm/d^ within the cloud; 
instead, the latter is a result of the magnetostatic 
calculation. Our method is to specify the shape 
of the cloud boundary Zci{r). On the boundary, 
the pressure Po is a known constant, while the po- 
tential, tpo = ^p[Zci{r)], is obtained through con- 
current solution of Poisson's equation. Prom the 



definition of q, equation (5), this quantity is there- 
fore determined everywhere within the cloud by 



g($) = Po exp(-V'o/a^), 



(16) 



where ipo is known along each fiux tube. With 
successive iterations, estimates of g, ip, and A im- 
prove, until convergence is achieved. (Por conver- 
gence criteria and other numerical details, sec the 
Appendix.) Henceforth, we refer to this technique 
as the q-method. We revert to the customary, or 
free-boundary method when constructing equilib- 
ria that have evolved from the initial state. 

2.4. Quasi-static Evolution 

At typical molecular densities of 10^ — 10^ cm~^, 
the cloud can remain in equilibrium with a frozen 
field for roughly 10^ yr. Over longer periods, 
however, self-gravity causes the neutral species to 
drift inward relative to the ions, which are tied to 
the magnetic field linos (Mestel & Spitzer 1956). 
The drift velocity of neutrals with respect to ions, 
V(j = V — Vj (no subscript indicating neutrals), de- 
pends both on the cloud's level of ionization and 
on the amount of collisional drag between neutrals 
and ions. At the low density contrasts character- 
izing our equilibria, |v| <C a, so that the p D'v/Dt 
term normally appearing on the right-hand side 
of the equation of motion (equation 1) can safely 
be ignored. Wc arc thus describing the quasi- 
static phase of evolution (Nakano 1984), during 
which the cloud slowly progresses along a sequence 
of exact, magnetostatic equilibria. As the cloud 
evolves, it becomes more centrally condensed, and 
gravity becomes relatively more important than 
magnetic forces near the cloud center. Once the 
neutral velocities there approach the sound speed, 
the quasi-static approximation breaks down, and 
the cloud undergoes essentially hydrodynamic col- 
lapse (e.g., Piedler & Mouschovias 1993). In this 
paper, we shall follow the quasi-static phase only, 
while also keeping track of the neutral velocities 
during the evolution. 

The drift speed is given by 

Vd = V- Vi = --(V X B) X B 
Pi 

= -(7C)-V-'/' Jl exp(-^/a2)V$, (17) 

where we have chosen the opposite sign convention 
to Nakano (1979), who gave an equivalent expres- 
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sion for the magnetic field (or ion) drift. Here 
Ti = (7/c)~^ is the damping time of ions relative 
to neutrals, pi is the mass density of the ions, and 
7 — 4.28 X 10"'^^ cm'' g~^ is the frictional drag 
coeflBcient (Nakano 1979). In the second step, we 
have made use of equations (3) and (7), and have 
assumed pi = Cp^l"^, where C = 4.46 x 10"-^^ g-^/^ 
cm-3/2 (Nakano 1979; Elmegreen 1979). This ex- 
pression for Pi assumes that ionization from cos- 
mic rays balances recombination. We ignore neg- 
atively charged grains, since their effect is small 
for the typical neutral gas densities we consider in 
this paper (n <C 10^ cm~^). For a poloidal mag- 
netic field that decreases outward in the cloud, 
both dqjd^ and d^/dr are positive, so Vd,r < 0, 
indicating inward drift of the neutrals (note that 
\vd,z\ {"fdiA during most of the quasi-static 
phase). The neutral velocity, v, is calculated us- 
ing the prescription of Lizano & Shu (1989).'* The 
qualitative behavior of both and v is shown in 
Figure 2. 

Neutral-ion drift leads to a redistribution of 
mass with magnetic flux in the cloud. More pre- 
cisely, the time rate of change of the mass con- 
tained within a given (axisymmetric) surface of 
constant magnetic flux <I> is 



= - P Vd-dS, 



(18) 



where the surface integration is performed along 
the entire flux tube $ = constant (Figure 2), and 
where the definition (14) of the mass function has 
been extended to the quasi-static case, i.e. m = 
m{^,t). Since dS points outward by definition, 
generally opposite to the direction of v^, dm/dt 
is usually positive for all $ < $o- It then follows 
that $0 is usually a decreasing function of time, 
i.e., the cloud as a whole loses flux. In the present 
coordinates, equations (17) and (18) combine to 
give 



d_ 
dt 



4n dq 
5$ 



I 



Zcl(*) 



"i/2exp(-V/a2) 



X r 



dr 



dz 



dz. 



(19) 



^Notc, however, that Lizano & Shu (1989) used values of 
7 and C which result in drift velocities 80 percent larger 
than those found in Tomisaka ct al. (1990) and the present 
paper. Their evolutionary timescales are correspondingly 
shorter. 



The quantity dm/dt vanishes at both = and 
$ = $0- At each time step At, one may adjust 
the mass-to-flux distribution as follows: 



dm 



{<^,t + At) = 



+ 



dm 



9$ 



dm , ^ , 



.(20) 



Given dm/d^ at the next time step, q{^,t + At) 
may be calculated from equation (15). Recall that 
{dm/ d^)t=o is known from the initial state, con- 
structed via the g— method (equation 13). During 
the subsequent evolution the boundary shape is al- 
lowed, to vary, with the solution at all times obey- 
ing the same boundary conditions as in the mag- 
netostatic case. Equation (15) is therefore used 
to calculate q{^) during the quasi-static phase of 
the calculation, and the new equilibrium is con- 
structed via the free-boundary method. Here, the 
boundary location Zc\{^) is calculated by finding 
the locus of points where /o($, z) = po = Pq/o'^. 

2.5. Nondimensionalization and Free Pa- 
rameters 

2.5.1. Magnetostatic equations 

We nondimensionalize our equations with re- 
spect to the sound speed a, cloud boundary den- 
sity Po, and background magnetic field i?oo- From 
these, we construct the following dimensionless 
quantities for the remaining variables: 

r' = r • (47rG/9o)^^V«> 
z' = z- {A-KGpof^/a, 

A' = A-{A^Gpof'^/aB^, 

$' = $-(47rGpo)/a'Soo, 
q' = q/a^po, 

M' = M- {AnGpof/'^/a^po, 
t' = t- {AT^Gpof^. 

Observations suggest that the typical cloud core 
temperature is T ^ 10 K (Jijina et al. 1999), 
while a typical number density outside cores is 
uh^ ~ 10^ cm~^ (Nercessian et al. 1988). Using 
a = {ksT / pY^'^ with p = 2.33 mH, one then finds 
that r' = 1 corresponds to a dimensional length of 
Lq = 0.11 pc. Similarly, M' = 1 corresponds to 
Mo = 0.070 Mq, and t' = 1 to io = 5.5 x 10^ yr. 
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With these definitions, the dimcnsionless forms 
of the fundamental equations (8) and (9) are 



d_ 
dr' 



r' dr' 
r 



+ 



r' dr' 



^ expl-V^ ) — 




dr' J ^ 9z'2 
q' exp(-^') 




(interior) (21) 
(exterior), 



(interior) (22) 
(exterior). 



The dimensionless parameter a appearing on the 
right-hand side of equation (21) is defined in terms 
of fiducial quantities by 

That is, a is the ratio of magnetic pressure far from 
the cloud to gas pressure at the cloud surface. 

2.5.2. Quasi-static evolution equations 

The dimensionless forms of equations (17) and 
(19) are 

v;, = ~C, ^ eM-i^') (V$)' (23) 



and 



dm' 



= C2 



dq' 



Zcl(*') 



exp(-V') 



, d^ 
'd7 



1 



dr' 
'd? 



dz', (24) 



where the numerical coefficients on the right- 
hand sides of equations (23) and (24) are de- 
fined in terms of dimensional quantities by Ci = 
{AnGf^hC = 0.0480 and C2 = 47rCi = 0.603, 
respectively. 

The point-mass boundary condition (equation 
11), in dimensionless form, becomes 



^' = 



-M' 



47r(r'2 -Hz'2)i/2 



at r' = R', z' = Z'. 



(25) 

The remaining boundary conditions are just the 
primed forms of those in §2.2, upon making the 
replacement = 1. 



2.5.3. Free parameters and functions 

To eonstruc;t ciquilibria under the point-mass 
boundary condition, we need to specify one pa- 
rameter, a, and one function, the cloud shape, 
Zc\{r). Measurements of magnetic field strengths 
in low-mass star-forming regions (mostly upper 
limits) suggest a mean of roughly 20 /xG (Crutcher 
1999), corresponding to a = 11. Given the un- 
certainty in the field measurements, we have per- 
formed calculations for a = 1 and a = 10. 

We constructed both spherical clouds and 
oblate/prolate ellipsoidal clouds with shapes given 
by 



/2\l/2 
) ^ 



(26) 



where R'q and Zq = Z'^^iO) are the principal axes of 
the ellipsoid. In the oblate and prolate cases, the 
cloud's axial ratio Z'q/R'q was fixed at 1/2 and 2, 
respectively. These figures are in accord with the 
most likely intrinsic mean value derived from ob- 
servations (Myers et al. 1991; Ryden 1996). As we 
describe below in §3, equilibrium sequences were 
constructed by considering a range of sizes R'q. 

Finally, under the periodic boundary condition, 
we need to specify the additional parameter 2Z , 
the intercore spacing. For the typical spacings in 
Taurus discussed in §1, i.e. 0.50 pc ^ 2Z ^ 1.0 pc, 
wc have 2.3 <^ Z' <^ 4.5. We adopt the value 
Z' = 2.9, corresponding to a dimensional spacing 
of 2Z = 0.65 pc. We also fix the extent of the 
computational volume in the r-direction at R' = 

3. INITIAL STATES 

In this section, we present the entire set of equi- 
libria obtained using the q'— method (§2.3). We re- 
fer to them as "initial" states, since we shall later 
be concerned with their quasi-static evolution via 
ambipolar diffusion. 

3.1. Calculation of Equilibrium Sequences 

For a given cloud shape and boundary condi- 
tion (point- mass or periodic), a sequence of states 
is obtained, each member of which is uniquely 
specified by its center-to-surface density contrast 
Pc- (Henceforth, we shall use only nondimensional 
quantities, omitting the primes, unless otherwise 
specified.) As pc increases, Rq and M first rise 
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to maximum values, i?o.max and A{f,„ax, respec- 
tively. Both quantities then decline. For much 
larger pc than covered in this study, Rq and M 
should undergo damped, oscillatory behavior, as 
found in the sequence of isothermal spheres (Chan- 
drasekhar 1939). 

For the low-/9c portion of each sequenc;e, we used 
the radius i?o as an input parameter for the mod- 
els. That is, we specified the cloud shape and 
then used the ij-method to find both pc and the 
internal structure of each model. At this stage, 
convergence was not overly sensitive to the initial 
guesses for the gravitational and magnetic fields. 
In practice, tpif^^) = and A{r,z) = r/2 proved 
satisfactory initial guesses, where the latter cor- 
responds to a uniform magnetic field set to the 
background value. 

In each sequence, the equatorial and polar ra- 
dius always peaked well before the cloud mass. 

This feature is familiar from the sequence of 
isothermal, nonmagnetic (Bonnor-Ebert) spheres, 
for which i?o,max = 1-822 and M^ax = Mbe = 
52.66 in our units. In general, one cannot pro- 
ceed to models with higher pc than that of i?o,max 
by specifying the cloud shape. This is because, 
along this portion of the sequence, each radius 
i?o < -Ro.max may correspond to more than one 
Pc— value and internal structure. We therefore 
modified the procedure, with the technique de- 
pending on the particular sequence (see below). 
Proceeding in the appropriate manner to larger 
Pc one eventually surpasses the mass peak of the 
sequence, Mmax- Since this peak signifies, at least 
roughly, the transition to dynamically unstable 
clouds (§3.4), we halted the search for equilibria 
soon thereafter. 

3.2. Isolated Clouds 

3.2.1. Force-Free States 

In the absence of a magnetic field and under 
the point-mass boundary condition, all equilibria 
are perfect spheres. The entire class of such iso- 
lated configurations bounded by a constant surface 
pressure was found by Ebert (1955) and Bonnor 
(1956). The addition of a magnetic field satisfying 
our prescribed boundary conditions (eq. 10) does 
not change this situation; i.e., each sphere is still 
an exact, equilibrium state. The reason is that 
a uniform, parallel field at the background value 



exerts no internal forces. Thus, isolated, spherical 
clouds should be magnetically force-free states. 

As a first test of the method, we used it 
to construct the spherical sequence. For Rq < 
-Ro.max, initial guesses of ■^(r, z) = 0, A{r, z) = r/2 
led quickly to convergence. For the higher— pc 
models beyond the radius peak, we again chose 
A{r,z) = r/2 and used the exact, nonmagnetic 
solutions with higher pc as initial guesses for 
tp. Members of the resulting sequence all con- 
tained uniform fields. For example, we found that 
Be = I B (0,0)1 never differed fractionally from 
unity (i.e., the background value) by more than 
2 X 10~^. We also checked that the density profile 
of each configuration matched that of its nonmag- 
netic, isothermal counterpart having the same cen- 
tral density pc. Thus, spherical clouds under the 
point-mass boundary condition are indeed mem- 
bers of the Bonnor-Ebert sequence. 

3.2.2. Oblate Equilibria 

More general, nonspherical, configurations with 
shapes given by equation (26) require magnetic 
support. To offset the weight of the extra mass in 
the equatorial region, an oblate cloud has a field 
that rises inward, i.e., dB^/dr < for all r and z. 
We used the g-method to construct a number of 
states with axial ratio Zq/Rq = 1/2. Our purpose 
was mainly to compare results with those of previ- 
ous authors who employed the point-mass bound- 
ary condition (Mouschovias 1976a, b; Tomisaka et 
al. 1988a,b). 

Figure 3a shows a representative equilibrium 
with a = 10. As expected, the field lines have 
a modest inward bending toward the polar axis. 
The horizontal arrows in the figure represent the 
drift velocities in this initial state, and show the in- 
stantaneous pattern of mass redistribution in the 
cloud. All of these velocities point inward. The 
peak drift speed f(j,max = 0.033 occurs at the cloud 
equator (i.e., r = Rq, z = 0), where the magnitude 
of the current also has a maximum. 

Figure 36 displays the mass-to-flux distribution 
dm/d^ as a function of the magnetic flux Be- 
cause the density rises siibstantially toward the 
center, dm/d^ peaks toward the axis ($ = 0). 
We note again that previous authors specifled this 
function ab initio. For example, Tomisaka et al. 
(1988a) used the dm/d^ from a uniform-density 
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sphere threaded by a constant magnetic field. This 
dm/d^ also rises inward, but lacks the sharp cen- 
tral peak seen in our model. The resulting cloud is 
oblate, but more "boxy" in appearance than ours. 
(See Fig. 2d of Tomisaka et al. 1988a.) 

In constructing the entire sequence of oblate 
equilibria, it is necessary to include states with 
higher central density than that at the radius peak 
-Ro.max (§3-1). To find these higher— pc states, we 
proceeded as follows. Taking a previously con- 
verged state with Rq < .Ro.max, we first increased 
dm/d^ by a small, fixed amount at every This 
altered dm / d^ was then used in the free-boundary 
method (§2.3) to converge a new state. The re- 
sulting equilibrium, which was only roughly ellip- 
soidal, had the same total flux $o as the original 
cloud, but a slightly larger mass and central den- 
sity. Now using this state as an initial guess, its 
dm/d^ was increased again, resulting in a cloud 
of still higher mass and central density. The ex- 
ercise was repeated several times, until pc signif- 
icantly exceeded that of the last state converged 
via the g— method. This high-density, roughly el- 
lipsoidal model then served as an initial guess for 
the construction, now via the g— method, of a truly 
ellipsoidal state lying well past the radius peak. 
In a typical case, the highest-density guesses con- 
structed in this manner sufficed to converge equi- 
libria that lie not only beyond -Ro,max, but beyond 
Mmax as well. 

Figure 4 displays the masses of our oblate 
clouds as a function of their central density, for 
a-values of 1 and 10. Note that in these and 
all subsequent equilibrium sequences, the total 
flux threading the cloud is not fixed a priori, but 
varies continuously along the sequence. The dot- 
ted curve shows, for comparison, the equivalent 
plot for the spherical, force- free states. It is evi- 
dent that the mass of any oblate configuration is 
significantly greater than that of a force-free state 
with the same central density. For example, the 
a = 10 sequence has a peak mass of M = 93.3, 
attained at pc = 10.6. This mass is 1.77 times 
the corresponding spherical value, which occurs at 
Pc = 14.0. Note that while the central density cor- 
responding to the peak mass is certainly smaller 
than in the spherical case, the broadness of the 
peak and finite sampling of the numerical calcula- 
tions prevent an accurate determination of Pc by 
the present method; our estimate may be in error 



by as much as 5-10 percent. 

The increased mass of oblate configurations 
has been noted by other authors (beginning with 
Mouschovias 1976b), and is readily understood. 
As Figure 3a indicates, the actual bending of the 
magnetic field is relatively slight in all our models. 
Hence, force balance in the vertical direction in- 
volves mainly gravity and thermal pressure. Sup- 
pose now that we envision building up an oblate 
cloud by adding mass equatorially to an initially 
spherical configuration. Then pc, which is largely 
determined by the weight of the central column, 
rises only slightly even for significant mass addi- 
tion. By the same token, it is clear that configura- 
tions even more oblate than ours will have corre- 
spondingly higher masses at any central density.^ 

One feature of the oblate models that may 
appear somewhat surprising is that their masses 
are insensitive to the background magnetic field 
strength, as parameterized by a. This is again 
a consequence of the relatively minor role played 
by the magnetic field in the vertical equilibrium. 
When the vertical magnetic force is small, the po- 
lar radius is set almost entirely by pc, as can be 
seen in the top panel of Figure 5. This property, 
coupled with the fact that the cloud shapes are 
constrained to be identical, means that two clouds 
of the same pc must have nearly identical masses, 
whatever their values of a. 

The above results imply that it is the shape 
of a cloud, not the strength of the magnetic field 
threading it, that determines the maximum equi- 
librium mass. However, the degree of field line 
bending in a particular state does depend on the 
background field, as well as on the cloud shape. 
We remarked earlier that the extra equatorial 
mass in oblate clouds creates an additional gravi- 
tational force in the radial direction. In response, 
the magnetic field bends inward, until the asso- 
ciated tension helps offset this force. The lower 
panel of Figure 5 quantifies the situation by plot- 
ting the central field value Be as a function of pc- 
For a = 1, Be rises substantially, reflecting the in- 
creased field line bending in higher-mass clouds. 
This rise is much less for a = 10, where the 

Any extra mechanical support beyond thermal pressure will 
allow higher cloud masses. For the mass increase due to 
rotation, see Figure 9 of Stahler (1983) and Figure 9 of 
Kiguchi et al. (1987). 
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stronger field resists the equatorial gravity more 
effectively. 

3.2.3. Prolate Equilibria 

We next consider clouds that arc prolate ellip- 
soids, with their long axes parallel to the back- 
ground magnetic field. Imagine creating such an 
object by shaving off equatorial mass from an ini- 
tially spherical configuration. By our previous rea- 
soning, this change would only marginally affect 
the central density. Thus, a prolate cloud is less 
massive than a spherical one with the same pc- In 
addition, the pressure gradient is higher in the ra- 
dial direction than along the pole. Since the equa- 
torial gravity is now insufficient to offset this extra 
force, the magnetic field bends outward, creating 
inward tension. The central field Be is reduced 
from the background value, and dB^/dr is pos- 
itive. In this sense, isolated, prolate clouds are 
magnetically confined. 

Using the g-method, we constructed sequences 
of prolate clouds with axial ratio Zq/Rq = 2. A 
representative state from the a = 1.5 sequence 
is shown in Figure 6a. The outward bowing of 
the field lines is evident, and is particularly strong 
near the center. Here the field attains its minimum 
value, less than 10 percent of the background. 
The figure also shows the drift velocities associ- 
ated with ambipolar diffusion. These now point 
outward, reflecting the new fleld curvature. The 
maximum drift speed, here equal to 0.3, occurs at 
the cloud equator. We were not able to obtain 
equilibria with lower values of a, most likely due 
to the very large gradient of B at the origin as 
Be approaches zero. Note that, because Ud.max is 
a significant fraction of the sound speed in these 
low-a states, our basic assumption of quasi-static 
equilibrium begins to fail. Wc shall revisit this is- 
sue in §4, when we attempt to follow explicitly the 
temporal evolution of prolate configurations. 

Figure 66 shows the mass-to-flux distribution 
in the cloud discussed above. Comparing with the 
analogous Figure 36 for an oblate state, we see 
first that the maximum flux $o is considerably 
smaller. This reduction is a consequence of the 
prolate cloud's smaller cross section in the equato- 
rial plane. We also see that dm/d^ climbs much 
more steeply toward the central axis, a result of 
both the increased column of gas near the pole 
and the reduced flux in that region. 



We display the cloud mass M as a function of 
Pc for the a = 1.5 and a = 10 sequences in Figure 
7. As in the oblate case, the curves are rather in- 
sensitive to a. However, the deficit in equatorial 
mass now implies that M is significantly less than 
the force- free value at every pc. Specifically, the 
mass is always below the Bonnor-Ebert value. The 
a = 10 curve has the usual broad peak, now cen- 
tered at M = 37.8, Pc = 14.9. This mass is 0.72 
times the corresponding maximum in the spher- 
ical sequence. For a = 1.5, however, the central 
field becomes so small that we were unable to con- 
struct models for pc beyond 15.5. Prolate clouds 
thus have a maximum mass set by the background 
field strength, if the latter is sufficiently low. 

The top panel in Figure 8 shows the polar radii 
of all our configurations, again as a function of 
Pc. The insensitivity with respect to a is appar- 
ent, and demonstrates once more the minor role of 
the magnetic field in vertical force balance. The 
polar radii are now uniformly greater than their 
force-free counterparts. Finally, the bottom panel 
of Figure 8 shows in detail how the central field 
diminishes as the field lines bow outward. This 
curvature, and hence the drop in Be, is relatively 
small for a = 10. On the other hand. Be rapidly 
falls to zero in the a = 1.5 case, where the curva- 
ture is much more pronounced. 

3.3. Tidally Stressed Clouds 

3.3.1. Force- Free States 

Under the periodic boundary condition, mag- 
netically force-free configurations are no longer 
spherical. Our stipulation that dip/dz vanish at 
z = ±Z accounts for the gravitational influence 
of external masses located above and below the 
cloud of interest. Their tidal force distorts the 
cloud vertically, i.e., in the direction of the uni- 
form magnetic field. The resulting equilibria are 
mildly prolate, although not perfectly ellipsoidal. 
At fixed Z and R (the dimensions of the bound- 
ing cylindrical surface), there is a unique sequence 
parametrized again by pc- The cloud axial ratio 
now varies along this sequence. 

Figure 9 shows isodensity contours for three 
equilibria in a representative sequence. Here, we 
have chosen Z = 2.9 and i? = 10, following the 
discussion of §2.5.3. The clouds are nearly spher- 
ical at low Pc, as shown in the first panel. They 
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become maximally distorted when the cloud's po- 
lar axis Zq peaks. The central panel of Figure 9 
shows this pivotal configuration, at which Zq has 
reached 2.39, Zo/i?o = 1-43, and pc = 3.9. At 
higher central density, the polar axis shrinks and 
the cloud becomes more spherical again. The third 
panel is an example of such a high-density state. 

If we view the cloud as one in a periodic chain, 
then the value of Z represents half of the inter- 
cloud spacing. It is hardly surprising, then, that 
the character of our force-free sequences is sensi- 
tive to Z. The choice Z = 2.9, as we have seen, 
results in a sequence for which the maximum Zq 
is close to Z itself. For larger Z, the clouds are 
less tidally distorted. For Z <^ 2.8 and sufficiently 
large pc, the upper and lower boundaries of the 
cloud penetrate the cylindrical boundary, i.e., the 
individual clouds in the chain partially coalesce. 
Note that our results are much less sensitive to 
the value of ii, at least beyond some minimum. 
For the Z = 2.9 case, for example, we found that 
raising R to 20 had a negligible effect on our mod- 
els. Interestingly, the masses of the clouds in the 
force-free, tidal sequence do not differ significantly 
from their counterparts of the same central den- 
sity in the Bonnor-Ebert sequence; the agreement 
is better than 1 percent throughout the range of 
pc examined. 

3.3.2. Spherical Equilibria 

The properties of tidally stressed clouds of var- 
ious shapes may be deduced from their relation to 
the force-free states. Thus, consider a true prolate 
ellipsoid with a higher aspect ratio Zq/Rq than 
the mildly elongated, force-free state of the same 
Pc- The ellipsoidal configuration experiences less 
gravitational force toward its axis, and so has an 
internal magnetic field that bows outward. Its 
structure will be qualitatively similar to that of 
the prolate clouds studied in §3.2.3. By analogy, a 
highly oblate cloud has field lines that bow inward, 
again like its isolated counterpart. 

An interesting difference arises, however, when 
we study spherical clouds. These are no longer 
force-free. In fact, their aspect ratio is less than 
that of the corresponding force-free states. Their 
field lines thus bend inward, and dBz/dr < 0. 
During ambipolar diffusion, the drift velocities 
point toward the symmetry axis, leading to states 
of higher central density. We are thus motivated 



to study these objects further. Historically, vari- 
ous authors have considered spherical, magnetized 
clouds as a convenient idealization (Mestel 1965; 
Strittmatter 1966). However, these studies ig- 
nored the detailed field topology, as did the more 
recent evolutionary calculations of Safier, McKee, 
& Stabler (1996) and Li (1998). 

To construct the spherical sequence we again 
utilized the g-method, now with the tidal bound- 
ary condition, equation (12). As always, the pro- 
cedure was straightforward up to the radius peak. 
In converging states beyond this point, we used 
Bonnor-Ebert solutions of similar pc to supply an 
initial guess for V- We guessed A by first con- 
structing a slightly oblate cloud of nearly the same 
Pc under the same tidal boundary condition. Such 
a state has an inwardly bent magnetic field, simi- 
lar to the tidally stressed sphere of interest. As we 
moved to higher-pc spherical states, the amount of 
field-line bending increased. Hence, we needed to 
use progressively more oblate clouds to supply the 
guess for A. 

Figure 10a shows the isodensity contours and 
field lines for a cloud from the a = 10 sequence. 
The degree of field line bending and enhancement 
of Be are only slight for this relatively high a- 
value. As expected, the drift velocities point in- 
ward, with the largest velocity occurring at the 
equator. Both the field topology and pattern of 
drift speeds resemble those for the isolated, oblate 
cloud of Figure 3, which has nearly the same cen- 
tral density. The mass-to-flux distribution in the 
sphere, shown in Figure 106, is also qualitatively 
similar. 

Figure 11 displays the mass of our spherical 
clouds as a function of central density, for a = 1 
and 10. The outer cylindrical boundary was again 
held fixed at ^ = 2.9, i?, = 10. The force-free se- 
quence is represented by the dotted curve. For a 
given Pc, the magnetically supported states are ev- 
idently more massive than both force-free, tidally 
stressed clouds and Bonnor-Ebert spheres. In the 
first case, the greater mass is a consequence of the 
objects' larger equatorial extent, while the tidal 
boundary condition is responsible in the second. 
More specifically, the vanishing potential gradi- 
ent at z = ztZ effectively weakens self-gravity 
everywhere, so that a higher equilibrium mass is 
required at a given pc than for the point-mass 
boundary condition. The maximum mass of 60.1 
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in the tidally stressed, spherical sequence exceeds 
the force-free and Bonnor-Ebert values by 13 per- 
cent. Note finally that the mass in our sequence 
peaks earlier (pc = 9.9) than in cither the force- 
free (pc = 15.6) or Bonnor-Ebert {pc = 14.1) case. 

The top panel of Figure 12 shows the variation 
in radius along the a = 10 sequence. Notice that 
the radii are intermediate between the force-free 
polar and equatorial values, shown by the dashed 
and dotted curve, respectively. The radius curve 
for the a = 1 sequence is nearly identical to the 
a = W result, and so is not displayed. Both 
curves lie slightly above the radii of Bonnor-Ebert 
spheres (compare Figure 5). In the lower panel 
of the figure, we display the variation in the cen- 
tral magnetic field value. The two sequences are 
now clearly distinguishable, with the a = 1 curve 
having greater field bending and thus a higher Be- 

3.4. Remarks on Dynamical Stability 

Figures 4, 7, and 11 all show the cloud mass 
reaching a maximum value and then declining as 
Pc increases. Such behavior is well known in other 
contexts, where the mass extrema demarcate a 
stability transition within the fundamental mode 
of oscillation, and thus the onset of dynamical in- 
stability (Tassoul 1978). In addition to its stellar 
applications, this "static method" for diagnosing 
instability has also been used for rotating clouds 
(Stabler 1983; Kiguchi et al. 1987). 

The situation is more subtle for magnetostatic 
configurations. The physical basis of the static 
method is that two equal-mass equilibria on ei- 
ther side of an extremum may be considered end- 
states in a normal mode of zero frcciuency. For 
this interpretation to hold, however, the two con- 
figurations must have identical distributions of all 
quantities that are conserved during such an os- 
cillation. For example, two rotating, axymmetric 
states must have the same variation of specific an- 
gular momentum, since the oscillation exerts no 
internal torques. 

In the case of magnetized clouds, we may as- 
sume that flux freezing holds during the dynami- 
cal oscillation (even one of zero frequency). Thus, 
application of the static method requires that the 
mass-to-flux distribution not change along the se- 
quence. This is the case for magnetostatic mod- 
els constructed via the free-boundary method, i.e.. 



with a specified functional form of dm/d^ (§2.3). 
In our models, on the other hand, dm/d^ varies 
continuously along the sequence, so the static 
method is not applicable. While clouds of suffi- 
ciently high central density are certainly unstable, 
the actual transition must be displaced somewhat 
from the extremum. 

Another consideration comes into play for the 
prolate states. We have stressed that these equilib- 
ria, especially those of low a and high pc, are mag- 
netically confined, with the field bowing outward 
in the central region (recall Fig. 6). In laboratory 
plasmas, such a field curvature tends to be dy- 
namically unstable. If one interchanges two adja- 
cent flux tubes, the energy of the system decreases. 
(For an elementary demonstration, see Section 8.3 
of Nicholson 1983.) A normal mode analysis shows 
that the unstable perturbations have short wave- 
lengths perpendicular to the magnetic field, and 
long wavelengths parallel to it (Freidberg 1987, p. 
267). Self-gravity, neglected in the plasma con- 
text, should not materially affect the properties of 
such a mode. 

4. TIME EVOLUTION 

The presence of a non-zero drift velocity in the 
states constructed in the previous section indicates 
that they will evolve in time via ambipolar diffu- 
sion. Thus we next consider the quasi-static evolu- 
tion of these clouds, as outlined in §2.4. We chose 
as initial configurations the oblate, prolate, and 
spherical clouds depicted in Figs. 3, 6, and 10, 
respectively. The calculations in this stage em- 
ployed the free-boundary method, rather than the 
g'— method (§2.3 and 2.4). As soon as the cloud 
was allowed to evolve, its shape immediately be- 
gan to depart from a perfect ellipsoid. In those 
cases where the central concentration increased 
with time, the calculation terminated when the 
large density gradient near the cloud center made 
convergence difficult on our uniformly spaced grid. 
This typically occurred when the maximum neu- 
tral velocity was 0.2 0.5 times the sound speed. 
Hence the termination roughly coincided with the 
breakdown of the quasi-static approximation and 
the beginning of dynamical collapse. 
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4.1. Oblate Initial States 

We begin with the oblate cloud shown in Figure 
3, which has an axial ratio oi Zq/ Rq = 1/2. By 
use of the evolution equation (24), this state was 
evolved in time steps of magnitude At = 0.05, un- 
til the central density pc reached 115, at f = tfin = 
11.23. Beyond this point, we could find no more 
equilibria. We determined the cndpoint accurately 
by reducing the time step to At = 0.01 after 
t = 10. The final time corresponds to 6.2 Myr, us- 
ing the fiducial parameter values of Section 2.5.1. 

Figure 13 shows the cloud at the final converged 
time, tfin- Note the pronounced dip in the cloud 
boundary near the pole, a feature seen in previ- 
ous oblate, magnetized models (see, e.g., Figure 1 
of Tomisaka et al. 1990). This is a result of the 
increasing self-gravity in the central region, where 
gravity is unopposed by magnetic forces along the 
field lines. The neutral velocity increases toward 
the center, and the vectors become more radial. 
The maximum value of |v| is 0.45, and occurs near 
the midplane at {r,z) = (0.41,0.07). 

The increase of the cloud's central density and 

mass-to-flux ratio with time is shown in Figure 14. 
As matter drifts into the central region, (dm/d^)c 
rises steadily, increasing by 73 percent in the 
course of the evolution. Simultaneously, the cen- 
tral field lines squeeze together, causing the value 
of Be to more than double by tfin. Since the total 
flux decreases, the global distribution of dm,/d^ is 
becoming both narrower and more sharply peaked. 
The rapid increase of Pc, along with the large value 
of |v| near the center, both indicate that the quasi- 
static stage of evolution is ending. 

Finally, consider the evolution of the den- 
sity proflles. At all times, the profiles exhibit 
a turnover, leading to an envelope with nearly 
power-law density. The radius of this turnover, 
i.e., of the inner core region, shrinks with time, 
while the profile in the envelope steepens. Both 
results are in accord with earlier studies (e.g., 
Mouschovias 1991). 

4.2. Prolate Initial States 

The unique feature of prolate structures is that 
they are magnetically confined. That is, the field 
lines bow outward from the central axis, creating 
the tension needed to oppose the larger pressure 
gradient in the radial direction. Since the drift 



velocities also point outward, we expect the evo- 
lution through ambipolar diffusion to differ quali- 
tatively from the oblate case. 

We attempted to evolve the prolate initial state 
with a = 1.5 shown in Figure 6a, but did not suc- 
ceed. As discussed in §3.2.3, this state possesses a 
steep gradient of B at the origin, which makes con- 
vergence difficult using the free-boundary method. 
Since this gradient is not as extreme for larger val- 
ues of a, we followed instead the evolution of a 2:1 
prolate cloud with a = 10. Figure 15 shows the 
results. The upper left panel displays the initial 
state. At early times, the outward neutral veloc- 
ity causes the equatorial region to expand. By the 
time t = 5 (2.8 Myr with our fiducial parameters), 
the central density has dropped to 11.2 and the 
maximum neutral velocity has decreased to less 
than half its initial value (see Fig. 156). During 
the same interval, the magnetic field has begun to 
straighten; its central value has increased slightly 
from 0.89 to 0.90. 

By the next time displayed in Figure 15c, t = 20 
(11 Myr), the equator has moved out still further, 
while the pole has shrunk noticeably from its ini- 
tial position. The cloud is clearly becoming less 
prolate, and now has an axial ratio of Zq/ Rq = 1.3. 
The trends noted above in pc, B^, and i^max con- 
tinue. Because the field is straightening, ambipo- 
lar diffusion begins to slow down. By t = 50 (28 
Myr; Fig. 15d), where we stopped the calcula- 
tion, the cloud is approaching a spherical shape 
{Zq/Rq = 1.07), albeit slowly, since the maximum 
neutral velocity is only 0.5 percent of the sound 
speed. The central density has decreased to one- 
quarter of its starting value, so that the density 
profile is now much flatter than the initial one (see 
below). Over most of the cloud, the neutral veloc- 
ities are oriented perpendicular to the field lines, 
and still point outward. Taken together, these 
developments indicate that the configuration will 
never reach the point of dynamical collapse. 

Instead, our calculation suggests that the end- 
state of a prolate, magnetized cloud is a sphere 
threaded by a uniform magnetic field. Since a uni- 
form field exerts no forces, the sphere would be 
indistinguishable from a non-magnetic, pressure- 
bounded configuration of the same mass as the 
initial prolate cloud. To test this idea, we con- 
structed a sphere having the same central density 
{Pc = 3.82) as the prolate final state. The cor- 
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responding mass of this non-magnetic state is (in 
our units) M = 41.0, which is only 5 percent larger 
than that of the state shown in Figure 15d. Given 
that the axial ratio of this "final" prolate state is 
still 7 percent greater than unity, this difference 
is understandable. Finally, we compare the mean 
of the equatorial and polar radii, [RpR^Y/^ , with 
the radius of an equilibrium sphere of the same 
Pc- We find that the two quantities differ by only 
2 percent. 

The evolution of the prolate cloud's central den- 
sity and mass-to-flux ratio are shown in Figure 
16. The central density decreases smoothly from 
its initial value, asymptotically approaching the 
spherical, force-free value of pc = 3.8 at large 
t. The central mass-to-flux ratio also decreases 
steadily with time. Once the final state is reached, 
{dm/d^)c has decreased to one-third of its original 
value. Comparison of Figures 14 and 16 shows how 
the oblate and prolate clouds evolve in essentially 
opposite directions. Figure 17 displays the equato- 
rial and polar density profiles of the initial, prolate 
state, and the same profiles in the final configura- 
tion. We also show the profile of a non-magnetic 
sphere with a central density of 3.8. In contrast to 
the oblate case, the core region expands with time, 
with the envelope gradient becoming significantly 
shallower by t = 50. 

4.3. Spherical Initial States 

An isolated spherical cloud, i.e., one in which 
the gravitational potential obeys the point-mass 
boundary condition, is actually a member of the 
Bonnor-Ebert sequence (§3.2.1). Since the inter- 
nal magnetic field is uniform and force-free, the 
drift velocity should be identically zero, and the 
cloud should not evolve under ambipolar diffu- 
sion. We checked that this was the case. Thus, 
a cloud with an initial pc of 14.0 increased its cen- 
tral density by only 2 percent over a period of 10 
Myr. Its central magnetic field changed by 0.3 
percent during this time, remaining very close to 
unity, while the maximum neutral velocity never 
exceeded 0.012. Finally, the equatorial and polar 
radii never differed from the Bonnor-Ebert value 
by more than 2 percent. 

The more interesting case is a spherical cloud 
constructed with the tidal boundary condition. 
Here, the magnetic field bows inward, leading to 
drift velocities that point toward the central axis. 



Starting with the configuration shown in Figure 
10, we found that the central density rose mono- 
tonically. We were able to follow the evolution 
until t = 20.15, corresponding to 11.2 Myr in our 
fiducial units. By this time, Pc had increased from 
its initial value of 11.7 to 111. Note that a compa- 
rable rise in pc for the isolated, oblate cloud took 
about half the time, as the drift velocities were 
roughly twice as large. 

Imposition of the tidal boundary condition ef- 
fectively weakens the gravitational force in the ver- 
tical direction. Nevertheless, the climb in central 
density causes the polar region to compress, just 
as in the isolated, oblate case. The pronounced 
polar dip may be seen in Figure 18, which shows 
the density contours and magnetic field lines at 
t = 20.15. The pattern of drift velocities is also 
similar to the oblate case. 

Figure 19 displays the evolution of the cloud's 

central density and mass-to-flux distribution. By 
the time the final state has been reached, pc has 
increased by almost an order of magnitude. Be 
has nearly doubled, and the central mass-to-flux 
value has increased by 74 percent. Once again, pc 
is undergoing a very sharp increase, a sign that 
the dynamical phase of evolution is close at hand. 
The largest neutral speed, which occurs near the 
cloud center at {r,z) ~ (0.33,0.07), has the value 
0.16. The density profiles of the evolving configu- 
ration exhibit the same qualitative features as in 
the oblate cloud evolution, i.e., a shrinking core 
region and a steepening envelope gradient. 

5. DISCUSSION 

In this paper, we have taken the novel ap- 
proach of specifying the cloud shape, rather than 
the mass-to-fiux distribution, in order to construct 
equilibrium models. One of our significant find- 
ings is that these structures can be either oblate 
or prolate. Furthermore, we did not encounter, in 
our admittedly limited parameter survey, any re- 
strictions on the aspect ratio in cither case. This 
result contrasts sharply with the traditional view, 
inherited from Mouschovias (1976b), that magne- 
tostatic equilibria are necessarily flattened along 
the direction of the ambient field. We now see 
that this conclusion arose from the method of con- 
structing models. Specifically, mass-to-flux distri- 
butions that plateau toward the cloud center yield 



14 



flattened equilibria, while more centrally-peaked 
distributions result in elongated clouds.^ 

We are not the first, however, to construct 
elongated equilibria. Tomisaka (1991) added a 
toroidal field component to the purely poloidal 
models of Tomisaka et al. (1988a, b), and thereby 
obtained moderately prolate configurations. In 
these models, the toroidal field exerts a radial 
pinch that aids gravity in squeezing gas toward 
the central axis. Fiege & Pudritz (2000) have 
recently explored a large number of such mod- 
els. Their clouds, which span a wide range of 
aspect ratios, are all of low mass and density con- 
trast (M ^ 0.8 Mbe, Pc ^ 11, where Mbe is 
the maximum mass of the force-free spherical se- 
quence). Interestingly, a subset of their equilib- 
ria exhibit bowed-out field lines similar to those 
in our own prolate models (sec, e.g., their Fig- 
ures 2 and 7). In their case, however, the bow- 
ing is most pronounced near the cloud boundary, 
and only occurs when the toroidal field strength is 
comparable to, or exceeds, the poloidal one. Be- 
cause these twisted-field solutions have a non-zero 
toroidal component at infinity, the elongation of 
the clouds depends in part on the magnetic field 
configuration throughout the larger, parent body. 
In our solutions, the prolate shape arises solely 
from the local mass-to-flux distribution. 

Our second major result is that the fate of 
a cloud under ambipolar diffusion is determined 
both by its density contrast and by its shape. Fig- 
ure 20 summarizes this finding for the case of iso- 
lated clouds. Here we have reproduced, from Fig- 
ures 4 and 7, the mass curves for oblate, prolate, 
and spherical (force-free) configurations, where 
the first two have a — 10. The three open circles 
represent initial states whose evolution we followed 
numerically. 

Consider first the two oblate configurations. 
The initial state of higher density contrast has a 
mass larger than Mbe {M = 1.3 Mbe; Pc = 2.95). 
Since its field lines bend inward, the cloud's cen- 



'Both our models and all previous ones assume that the 
cloud is axisymmetric about the background field direction. 
Ward-Thompson et al. (2000) have recently published the 
first observations of magnetic fields in starless dense cores, 
using submillimeter polarimetry. In their three examples, 
the projected field direction appears to be at an oblique 
angle with respect to the cores' long axes. This result, if 
confirmed by future studies, may indicate that the assump- 
tion of axisymmetry requires modification (Basu 2000) . 



tral density increases as it evolves, and the repre- 
sentative point moves toward the right. Note that 
the contraction rate for this cloud is much slower 
than for the example in §4.1, which was at the 
mass peak of the oblate sequence. Initially, the 
equatorial region shrinks rapidly enough that the 
aspect ratio Zq/Rq actually increases, from 0.50 
at t = to 0.74 at t = 80. Eventually, however, 
the polar dip becomes so pronounced that Zq / Rq 
falls, reaching 0.55 at our last computed model 
{t = 119). At this time, the central density is 
Pc = 110 and is climbing rapidly. Thus, all oblate 
equilibria ultimately flatten before they collapse. 

By contrast, the initially oblate cloud of lower 
density contrast shown in Figure 20 has a mass be- 
low Mbe (M = 0.86 Mbe, Pc = 2.00). Although 
Pc still increases with time, the configuration does 
not collapse, but reaches the force-free, spherical 
state marked by the filled circle. By t = 100, the 
central density is 3.58, within 25 percent of the 
predicted final value of 4.80, and the aspect ratio 
is Zq/Rq = 0.94, very close to 1. Thus the cloud 
becomes less oblate and more spherical, while the 
field lines straighten out. 

Finally, we turn to the prolate sequence, also 
shown in Figure 20. Since the mass curve now lies 
wholly below the spherical one, no initial state will 
lead to dynamical collapse. The example shown 
by an open circle corresponds exactly to the state 
discussed previously in §4.2. There it was demon- 
strated that the central density falls as the con- 
figuration becomes more spherical, and that the 
cloud does not collapse. To summarize, all states 
move toward the force-free curve in the diagram, 
but only those with M < Mbe ^an actually reach 
the curve and so avoid collapse. An analogous 
situation holds for our spherical, tidally stressed 
equilibria. Referring again to Figure 11, we see 
that the initial state with pc — 11.7 lies above the 
force-free curve, and so inevitably collapses. 

Our third, and last, major result is that none 
of our initial equilibria, whether isolated or tidally 
stressed, evolves in a manner consistent with ob- 
served, star-forming clouds. As we reviewed in 
Section 1, the radio maps of dense cores are best 
explained as an ensemble of elongated, three- 
dimensional objects projected onto the plane of 
the sky. There is no appreciable difference between 
the shapes of cores with or without internal stars 
(Jijina et al. 1999). Thus, a correct theoretical ac- 
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count should start with an elongated configuration 
and have it remain elongated through the onset of 
dynamical collapse. 

The fault may lie in our neglect of external 
matter. In this study, we have followed the con- 
ventional route of surrounding the cloud by a ficti- 
tious, zero-density medium of finite pressure. Real 
dense cores, however, are embedded within larger 
filamentary structures. This fact alone hints why 
the initial states are elongated. During ambipolar 
diffusion, matter from the external reservoir can 
flow down magnetic field lines, piling up until it 
is halted by the internal pressure gradient. Thus, 
the polar flattening we obtained may be alleviated 
or even avoided entirely. It will be interesting to 
see how cloud evolution proceeds under these very 
different conditions. 
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A. APPENDIX: NUMERICAL TECHNIQUES 



To implement the rj— method, equations (21) and (22) were solved simultaneously on a imiform grid 
covering the upper right quadrant of the large cylinder; i.e., < r < R, < z < Z. We begin with an initial 
guess A^^\ '<P^^\ and a fixed cloud boundary Zc\{r). The function g($) is then given by equation (16), which 
we rewrite as 

(Al) 



9(o)($)=Pocxp[-V.(°V«^ 



with the superscript indicating the zeroth iteration. Recall that the subscript zero on tp and P denotes their 
boundary values. As mentioned in Section 2, the specification of q{(^) is sufficient to calculate the right-hand 
side source terms of equations (21) and (22). At each subsequent iteration, solutions for A and ij) subject 
to the boundary conditions (10) and (11) (point-mass) or (12) (tidal) were obtained to within an accuracy 
of 10~^. These were used as provisional iterates A^}\ i/'*^^ in generating subsequent guesses, as explained 
below. At each step, the updated ij) was then used in equation (Al) to generate a new q. 

The number of grid points in the radial and vertical directions, Nr and Nz, for the different cloud shapes 
adopted in Section 3 were: 

Ohlate: Nr = 81, Nz = 41 
Prolate: Nr = 41, Nz = 81 
Spherical: Nr = 61, Nz = 41. 

The number of field lines, Nb, was chosen equal to Nr. To ensure accurate derivatives near the cloud 

center, the flux value of the i^^ field line was chosen as $, = rf/2, i = 1,2, .... Nb- In practice, it was 
necessary to switch back and forth between two interwoven meshes: one fixed spatially (r, z), and another, 
defined by the current field line positions ($, z) (Mouschovias 1976a). 
The sequence of iterates was chosen according to the scheme: 

= A^^^ + (1 - Oa) 
and = 0^ ^("^ + {l-e^) ^1"+'^ , 

where Oa and 9^ are constant relaxation parameters, taken to be between 0.5 and 1 in most of our calculations. 
A solution was deemed acceptable if the conditions 



An+1) 



< e 



and 



(n+l) 



< e 



(A2) 



were satisfied simultaneously, with e = 5 x 10 ^. 

We typically found solutions within 50 iterations or less, using an initial guess ip^^\r, z) = 0, A^°\r, z) = 
r/2. We checked that the solutions were relatively insensitive to: (i) increasing the number of grid points 
in either direction; (ii) increasing the values of R and Z; and (iii) allowing the program to continue for up 
to twice the number of iterations required to satisfy conditions (A2). None of these resulted in a change of 
more than 5 percent in ijj or A at any grid point, or of more than 8 percent in the central density. 

The numerical implementation of the free-boundary method was discussed by Mouschovias (1976a). In 
applying the method, we imposed an additional convergence criterion on the location of the cloud boundary: 



An) 



< ez 



(A3) 



where = Az/2. Thus, at each timestep, convergence was achieved when conditions (A2) and (A3) were 
simultaneously satisfied. 

In order to evolve the cloud over the hundreds of timesteps typical of our calculations in Section 4, 
it was essential that the mass-to-flux distribution given by equation (20) remain smooth throughout the 
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evolution. To ensure this was the case, we approximated the function at each timestep by a least- 

squares polynomial. Thus, dq/d^ was also a polynomial, lending an additional degree of smoothness to 
equations (19) and (20). 

REFERENCES 

Basu, S. 2000, ApJ, 540, L103 

Benson, P.J., & Myers, P.C. 1989, ApJS, 71, 89 

Bonnor, W.B. 1956, MNRAS, 116, 350 

Chandrasekhar, S. 1939, An Introduction to the Study of Stellar Structure (New York: Dover) 

Crutcher, R.M. 1999, ApJ, 520, 706 

Dungey, J. 1953, MNRAS, 113, 180 

Ebert, R. 1955, ZAp, 37, 216 

Elmegreen, B.G. 1979, ApJ, 232, 729 

Fiedler, R.A., & Mouschovias, T.Ch. 1992, ApJ, 391, 199 

Fiedler, R.A., & Mouschovias, T.Ch. 1993, ApJ, 415, 680 

Ficge, J.D., & Pudritz, R.E. 2000, ApJ, 534, 291 

Freidberg, J. P. 1987, Ideal Magnetohydrodynamics (New York: Plenum) 

Harju, J., Walmsley, CM., & Wouterloot, J.G.A. 1993, A&AS, 98, 51 

Jijina, J., Myers, P.C, & Adams, F.C 1999, ApJS, 125, 161 

Kiguchi, M., Narita, S., Miyama, S.M., & Hayashi, C. 1987, ApJ, 317, 830 

Li, Z.-Y. 1998, ApJ, 493, 230 

Lizano, S., & Shu, F.H. 1989, ApJ, 342, 834 

McKee, CF., Zweibel, E.G., Goodman, A.A., & Heiles, C. 1993, in Protostars and Planets III, ed. E.H. 
Levy & J.I. Lunine (Tucson: University of Arizona Press), 327 

Mestel, L. 1965, QJRAS, 6, 265 

Mestel, L., & Spitzer, L. 1956, MNRAS, 116, 503 

Mouschovias, T.Ch. 1976a, ApJ, 206, 753 

Mouschovias, T.Ch. 19766, ApJ, 207, 141 

Mouschovias, T.Ch. 1991, in The Physics of Star Formation and Early Stellar Evolution, ed. C.J. Lada & 
N.D. Kylafis (Dordrecht: Kluwer), 449 

Myers, P.C. 1999, in The Origin of Stars and Planetary Systems, ed. CJ. Lada & N.D. Kylafis (Dordrecht: 
Kluwer), 67 

Myers, P.C, & Benson, P.J. 1983, ApJ, 266, 309 
Myers, P.C, & Goodman, A. A. 1988, ApJ, 326, L27 



18 



Myers, P.C., Fuller, G.A., Goodman, A.A., & Benson, P.J. 1991, ApJ, 376, 561 

Nakano, T. 1979, PASJ, 31, 697 

Nakano, T. 1984, Fund. Cosm. Phys. 9, 139 

Nercessian, E., Castets, A., Benayoun, J.J., Cernicharo, J. 1988, A&A, 189, 207 
Nicholson, D.R. 1983, Introduction to Plasma Theory (New York: Wiley) 

Onishi T., et al. 1996, ApJ 465, 815 
Ryden, B.S. 1996, ApJ 471, 822 

Safier, P., McKee, C.F., & Stahler, S. 1997, ApJ, 485, 660 
Stabler, S.W. 1983, ApJ, 268, 165 
Strittmatter, P.A. 1966, MNRAS, 132, 359 

Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton: Princeton University Press) 
Tomisaka, K. 1991, ApJ, 376, 190 

Tomisaka K., Ikeuchi, S., & Nakamura, T. 1988a, ApJ 326, 208 
Tomisaka, K., Ikeuchi, S., & Nakamura, T. 19886, ApJ, 335, 239 
Tomisaka, K., Ikeuchi, S., & Nakamura, T. 1990, ApJ, 362, 202 

Ward-Thompson, D., Kirk, J.M., Crutcher, R.M., Greaves, J.S., Holland, W.S., & Andre, P. 2000, ApJ, 537, 
L135 



This 2-column preprint was prepared with the AAS lAl^jX macros v5.0. 



19 



Z 




R r 

Fig. 1. — The computational volume is the upper right quadrant of a right circular cylinder of radius R, 
half-height Z. The cloud is an axisymmetric structure with boundary Zci(r), centered on the origin, and 
embedded in a region of zero density but finite pressure. Magnetic field lines in this figure are indicated by 
solid lines. 
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Fig. 2. — The physics of ambipolar diffusion. A portion of the cloud, of mass Am, is enclosed between two 
flux tubes, as shown. Over time, the neutral particles drift inward relative to the field with a velocity v^. 
This drift velocity differs from v, that of the neutrals in an inertial reference frame. 
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Fig. 3. — (a) An oblate initial state with a = 10. Solid curves are isodensity contours and dotted curves 
are magnetic field lines. The heavy solid curve is the cloud boundary, which has an axial ratio Zq/Rq 
of 1/2. Arrows indicate drift velocities, Vd = v — v^. The dimensions of the computational volume are 
R = 5.52, Z = 2.76. (b) Mass-to-flux distribution as a function of magnetic flux $ for the same cloud. 
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Fig. 4. Mass as a function of density contrast (solid curves) for isolated, 2:1 oblate equilibria. Results for 
both values of a are displayed. The corresponding force-free sequence of spherical (Bonnor-Ebert) equilibria 
is also shown (dotted curve). 
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Fig. 5. — (top) Polar radius versus density contrast for the oblate 2:1 sequence (solid curves). The corre- 
sponding force- free sequence is also shown (dotted curve), (bottom) Central magnetic field versus density 
contrast, along the same sequences. 
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Fig. 6. — (a) An isolated, 2:1 prolate equilibrium cloud with a = 1.5. The magnetic field strength at the 
origin, Be, is about 10 percent of its asymptotic value. Here, the drift speeds point outward at all points 

in the cloud. The maximum drift speed, Ud.max, is 30 percent of the sound speed. The dimensions of the 
computational volume are R = 2.24, Z = 4.4:7. (b) The mass-to-flux function of the cloud in (a). 



25 




Fig. 7. — Same as Figure 4, but for the sequence of 2:1 prolate equilibria (solid curve). The dotted curve 
shows the sequence of force-free (spherical) equilibria. 
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Fig. 8. — (top) Polar radius versus density contrast for the prolate 2:1 sequence (solid curves). The corre- 
sponding force- free sequence is also shown (dotted curve), (bottom) Central magnetic field versus density 
contrast, along the same sequences. 
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Fig. 9. Force- free equilibria under the tidal boundary condition. The dimensions of the computational 
volume are R = 10.0, Z = 2.9. (a) A low-density cloud with pc = 1.5. (b) Maximally distorted cloud with 
pc = 3.9. (c) A high-density cloud with pc = 100. 
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Fig. 10. — (a) A tidally stressed, spherical equilibrium cloud with a = 10. The dimensions of the computa- 
tional volume are R = 10.0, Z = 2.9. (b) The mass-to-flux function of the cloud in (a). 
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Fig. 11. — Same as Figure 4, but for the sequence of spherical equihbria under the tidal boundary condition 
(solid curve). The dashed curve shows the sequence of force- free (prolate) equilibria. 
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Fig. 12. — (top) Radius versus density contrast for spheres under the tidal boundary condition with a = 10 
(solid curves). The results are nearly identical for a = 1. The equatorial and polar radii of the force- free 
sequence are also shown, (bottom) Central magnetic field versus density contrast, along the a = 1 and 
a = 10 sequences. 
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Fig. 13. — Final converged state of the isolated, 2:1 oblate cloud shown in Fig. 3. Arrows indicate neutral 
velocities, v. The dashed curve shows the initial cloud boundary. 
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Fig. 14. — Time evolution of the central density (solid curve) and central mass-to-flux (dashed curve) for 
the oblate initial state of Fig. 3. 



33 




Fig. 15. — Time evolution of an isolated, 2:1 prolate cloud with a = 10. Panel (a) shows the density contours 
(solid curves) and magnetic field lines (dotted curves) in the initial state. Panels (b), (c), and (d) show the 
cloud at the times t = 5,20, and 50, respectively. Arrows indicate drift velocities in panel (a), neutral 
velocities in panels (b), (c), and (d). 
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Fig. 16. — Time evolution of the central density (solid curve) and central mass-to-flux (dashed curve) 
the initially prolate state of Fig. 15. 
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Fig. 17. — Density profiles of the prolate cloud in its initial and final states. Solid curves show the density 
along the midplane {z = 0); dashed curves, along the pole (r = 0). The density profile of an isolated, 
spherical cloud with the same central density is shown by a dotted curve. 
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Fig. 19. — Time evolution of the central density (solid curve) and central mass-to-flux (dashed curve) 
the spherical cloud of Fig. 10. 
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Fig. 20. — Mass as a function of density contrast for isolated equilibria considered in this paper (solid curves). 
Results are displayed for a = 10 only. Initial states are represented by open circles, final states by filled 
circles. 
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